Harvest index is a key trait for screening drought-tolerant potato genotypes (Solanum tuberosum)

Journal of Crop Science and Biotechnology

1 Setup

source("https://inkaverse.com/setup.r")
library(gtsummary)
library(factoextra)
library(corrplot)
library(emmeans)
library(multcomp)

devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.0 (2023-04-21 ucrt)
 os       Windows 11 x64 (build 22621)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  Spanish_Latin America.utf8
 ctype    Spanish_Latin America.utf8
 tz       America/Lima
 date     2023-05-19
 pandoc   2.19.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package       * version  date (UTC) lib source
 agricolae       1.3-5    2021-06-06 [1] CRAN (R 4.3.0)
 AlgDesign       1.2.1    2022-05-25 [1] CRAN (R 4.3.0)
 askpass         1.1      2019-01-13 [1] CRAN (R 4.3.0)
 boot            1.3-28.1 2022-11-22 [2] CRAN (R 4.3.0)
 broom.helpers   1.13.0   2023-03-28 [1] CRAN (R 4.3.0)
 cachem          1.0.8    2023-05-01 [1] CRAN (R 4.3.0)
 callr           3.7.3    2022-11-02 [1] CRAN (R 4.3.0)
 cellranger      1.1.0    2016-07-27 [1] CRAN (R 4.3.0)
 cli             3.6.1    2023-03-23 [1] CRAN (R 4.3.0)
 cluster         2.1.4    2022-08-22 [2] CRAN (R 4.3.0)
 codetools       0.2-19   2023-02-01 [2] CRAN (R 4.3.0)
 colorspace      2.1-0    2023-01-23 [1] CRAN (R 4.3.0)
 combinat        0.0-8    2012-10-29 [1] CRAN (R 4.3.0)
 corrplot      * 0.92     2021-11-18 [1] CRAN (R 4.3.0)
 cowplot       * 1.1.1    2020-12-30 [1] CRAN (R 4.3.0)
 crayon          1.5.2    2022-09-29 [1] CRAN (R 4.3.0)
 curl            5.0.0    2023-01-12 [1] CRAN (R 4.3.0)
 devtools      * 2.4.5    2022-10-11 [1] CRAN (R 4.3.0)
 digest          0.6.31   2022-12-11 [1] CRAN (R 4.3.0)
 dplyr         * 1.1.2    2023-04-20 [1] CRAN (R 4.3.0)
 DT              0.27     2023-01-17 [1] CRAN (R 4.3.0)
 ellipsis        0.3.2    2021-04-29 [1] CRAN (R 4.3.0)
 emmeans       * 1.8.6    2023-05-11 [1] CRAN (R 4.3.0)
 estimability    1.4.1    2022-08-05 [1] CRAN (R 4.3.0)
 evaluate        0.21     2023-05-05 [1] CRAN (R 4.3.0)
 factoextra    * 1.0.7    2020-04-01 [1] CRAN (R 4.3.0)
 FactoMineR    * 2.8      2023-03-27 [1] CRAN (R 4.3.0)
 fansi           1.0.4    2023-01-22 [1] CRAN (R 4.3.0)
 fastmap         1.1.1    2023-02-24 [1] CRAN (R 4.3.0)
 flashClust      1.01-2   2012-08-21 [1] CRAN (R 4.3.0)
 forcats       * 1.0.0    2023-01-29 [1] CRAN (R 4.3.0)
 fs              1.6.2    2023-04-25 [1] CRAN (R 4.3.0)
 gargle          1.4.0    2023-04-15 [1] CRAN (R 4.3.0)
 generics        0.1.3    2022-07-05 [1] CRAN (R 4.3.0)
 ggplot2       * 3.4.2    2023-04-03 [1] CRAN (R 4.3.0)
 ggrepel         0.9.3    2023-02-03 [1] CRAN (R 4.3.0)
 glue            1.6.2    2022-02-24 [1] CRAN (R 4.3.0)
 googledrive   * 2.1.0    2023-03-22 [1] CRAN (R 4.3.0)
 googlesheets4 * 1.1.0    2023-03-23 [1] CRAN (R 4.3.0)
 gsheet        * 0.4.5    2020-04-07 [1] CRAN (R 4.3.0)
 gt              0.9.0    2023-03-31 [1] CRAN (R 4.3.0)
 gtable          0.3.3    2023-03-21 [1] CRAN (R 4.3.0)
 gtsummary     * 1.7.1    2023-04-27 [1] CRAN (R 4.3.0)
 haven           2.5.2    2023-02-28 [1] CRAN (R 4.3.0)
 highr           0.10     2022-12-22 [1] CRAN (R 4.3.0)
 hms             1.1.3    2023-03-21 [1] CRAN (R 4.3.0)
 htmltools       0.5.5    2023-03-23 [1] CRAN (R 4.3.0)
 htmlwidgets     1.6.2    2023-03-17 [1] CRAN (R 4.3.0)
 httpuv          1.6.11   2023-05-11 [1] CRAN (R 4.3.0)
 httr            1.4.6    2023-05-08 [1] CRAN (R 4.3.0)
 huito         * 0.2.2    2023-01-24 [1] CRAN (R 4.3.0)
 inti          * 0.6.0    2023-01-24 [1] CRAN (R 4.3.0)
 jsonlite        1.8.4    2022-12-06 [1] CRAN (R 4.3.0)
 klaR            1.7-2    2023-03-17 [1] CRAN (R 4.3.0)
 knitr         * 1.42     2023-01-25 [1] CRAN (R 4.3.0)
 labelled        2.11.0   2023-04-11 [1] CRAN (R 4.3.0)
 later           1.3.1    2023-05-02 [1] CRAN (R 4.3.0)
 lattice         0.21-8   2023-04-05 [2] CRAN (R 4.3.0)
 leaps           3.1      2020-01-16 [1] CRAN (R 4.3.0)
 lifecycle       1.0.3    2022-10-07 [1] CRAN (R 4.3.0)
 lme4            1.1-33   2023-04-25 [1] CRAN (R 4.3.0)
 lubridate     * 1.9.2    2023-02-10 [1] CRAN (R 4.3.0)
 magick        * 2.7.4    2023-03-09 [1] CRAN (R 4.3.0)
 magrittr        2.0.3    2022-03-30 [1] CRAN (R 4.3.0)
 MASS          * 7.3-58.4 2023-03-07 [2] CRAN (R 4.3.0)
 Matrix          1.5-4    2023-04-04 [2] CRAN (R 4.3.0)
 memoise         2.0.1    2021-11-26 [1] CRAN (R 4.3.0)
 mime            0.12     2021-09-28 [1] CRAN (R 4.3.0)
 miniUI          0.1.1.1  2018-05-18 [1] CRAN (R 4.3.0)
 minqa           1.2.5    2022-10-19 [1] CRAN (R 4.3.0)
 mnormt          2.1.1    2022-09-26 [1] CRAN (R 4.3.0)
 multcomp      * 1.4-23   2023-03-09 [1] CRAN (R 4.3.0)
 multcompView    0.1-9    2023-04-09 [1] CRAN (R 4.3.0)
 munsell         0.5.0    2018-06-12 [1] CRAN (R 4.3.0)
 mvtnorm       * 1.1-3    2021-10-08 [1] CRAN (R 4.3.0)
 nlme            3.1-162  2023-01-31 [2] CRAN (R 4.3.0)
 nloptr          2.0.3    2022-05-26 [1] CRAN (R 4.3.0)
 openssl         2.0.6    2023-03-09 [1] CRAN (R 4.3.0)
 pillar          1.9.0    2023-03-22 [1] CRAN (R 4.3.0)
 pkgbuild        1.4.0    2022-11-27 [1] CRAN (R 4.3.0)
 pkgconfig       2.0.3    2019-09-22 [1] CRAN (R 4.3.0)
 pkgload         1.3.2    2022-11-16 [1] CRAN (R 4.3.0)
 prettyunits     1.1.1    2020-01-24 [1] CRAN (R 4.3.0)
 processx        3.8.1    2023-04-18 [1] CRAN (R 4.3.0)
 profvis         0.3.8    2023-05-02 [1] CRAN (R 4.3.0)
 promises        1.2.0.1  2021-02-11 [1] CRAN (R 4.3.0)
 ps              1.7.5    2023-04-18 [1] CRAN (R 4.3.0)
 psych         * 2.3.3    2023-03-18 [1] CRAN (R 4.3.0)
 purrr         * 1.0.1    2023-01-10 [1] CRAN (R 4.3.0)
 questionr       0.7.8    2023-01-31 [1] CRAN (R 4.3.0)
 R6              2.5.1    2021-08-19 [1] CRAN (R 4.3.0)
 rappdirs        0.3.3    2021-01-31 [1] CRAN (R 4.3.0)
 Rcpp            1.0.10   2023-01-22 [1] CRAN (R 4.3.0)
 readr         * 2.1.4    2023-02-10 [1] CRAN (R 4.3.0)
 remotes         2.4.2    2021-11-30 [1] CRAN (R 4.3.0)
 rlang           1.1.1    2023-04-28 [1] CRAN (R 4.3.0)
 rmarkdown       2.21     2023-03-26 [1] CRAN (R 4.3.0)
 rstudioapi      0.14     2022-08-22 [1] CRAN (R 4.3.0)
 sandwich        3.0-2    2022-06-15 [1] CRAN (R 4.3.0)
 scales          1.2.1    2022-08-20 [1] CRAN (R 4.3.0)
 scatterplot3d   0.3-44   2023-05-05 [1] CRAN (R 4.3.0)
 sessioninfo     1.2.2    2021-12-06 [1] CRAN (R 4.3.0)
 shiny         * 1.7.4    2022-12-15 [1] CRAN (R 4.3.0)
 showtext        0.9-6    2023-05-03 [1] CRAN (R 4.3.0)
 showtextdb      3.0      2020-06-04 [1] CRAN (R 4.3.0)
 stringi         1.7.12   2023-01-11 [1] CRAN (R 4.3.0)
 stringr       * 1.5.0    2022-12-02 [1] CRAN (R 4.3.0)
 survival      * 3.5-5    2023-03-12 [2] CRAN (R 4.3.0)
 sysfonts        0.8.8    2022-03-13 [1] CRAN (R 4.3.0)
 TH.data       * 1.1-2    2023-04-17 [1] CRAN (R 4.3.0)
 tibble        * 3.2.1    2023-03-20 [1] CRAN (R 4.3.0)
 tidyr         * 1.3.0    2023-01-24 [1] CRAN (R 4.3.0)
 tidyselect      1.2.0    2022-10-10 [1] CRAN (R 4.3.0)
 tidyverse     * 2.0.0    2023-02-22 [1] CRAN (R 4.3.0)
 timechange      0.2.0    2023-01-11 [1] CRAN (R 4.3.0)
 tzdb            0.4.0    2023-05-12 [1] CRAN (R 4.3.0)
 urlchecker      1.0.1    2021-11-30 [1] CRAN (R 4.3.0)
 usethis       * 2.1.6    2022-05-25 [1] CRAN (R 4.3.0)
 utf8            1.2.3    2023-01-31 [1] CRAN (R 4.3.0)
 vctrs           0.6.2    2023-04-19 [1] CRAN (R 4.3.0)
 withr           2.5.0    2022-03-03 [1] CRAN (R 4.3.0)
 xfun            0.39     2023-04-20 [1] CRAN (R 4.3.0)
 xml2            1.3.4    2023-04-27 [1] CRAN (R 4.3.0)
 xtable          1.8-4    2019-04-21 [1] CRAN (R 4.3.0)
 yaml            2.3.7    2023-01-23 [1] CRAN (R 4.3.0)
 zoo             1.8-12   2023-04-13 [1] CRAN (R 4.3.0)

 [1] C:/Users/floza/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.0/library

──────────────────────────────────────────────────────────────────────────────

2 Googlesheets connect

url <- "https://docs.google.com/spreadsheets/d/1dfgpmCKdPmxRHozrZp0iE_xMGsvKTIcztDpMWYSEGaY"
# browseURL(url)
gs <- as_sheets_id(url)

3 Trait rename

geno <- c("CIP720088" = "G01"
          , "CIP392797.22" = "G02"
          , "CIP397077.16" = "G03"
          , "CIP398192.213" = "G04"
          , "CIP398180.612" = "G05"
          , "CIP398208.704" = "G06"
          , "CIP398098.119" = "G07"
          , "CIP398190.89" = "G08"
          , "CIP398192.592" = "G09"
          , "CIP398201.510" = "G10"
          , "CIP398203.244" = "G11"
          , "CIP398203.5" = "G12"
          , "CIP398208.219" = "G13"
          , "CIP398208.33" = "G14"
          , "CIP398208.620" = "G15")

traits <- c("SPAD_29" = "Chlorophyll concentration (SPAD) at 29 dap"
            , "SPAD_59" = "Chlorophyll concentration (SPAD) at 59 dap"
            , "SPAD_76" = "Chlorophyll concentration (SPAD) at 76 dap"
            , "SPAD_83" = "Chlorophyll concentration (SPAD) at 83 dap"
            , "HGT" = "Plant height (cm)"
            , "RWC" = "Relative water content (%)"
            , "LOP" = "Leaf osmotic potential (MPa)"
            , "LDW" = "Leaf dry weight (g)"
            , "SDW" = "Stem dry weight (g)"
            , "RDW" = "Root dry weight (g)"
            , "TDW" = "Tuber dry weight (g)"
            , "NTUB" = "Tuber number (N°)" 
            , "TRS" = "Total transpiration (mL)"
            , "LFA" = "Leaf area (cm^2^)"
            , "RTL" = "Root length (cm)"
            , "TDB" = "Total dry biomass (g)"
            , "HI"  = "Harvest index (HI)" 
            , "SLA" = "Specific leaf area (cm^2^ g^-1^)"
            , "RCC" = "Relative chlorophyll content (RCC)"
            , "WUE_B" = "Biomass water use efficiency (g L^-1^)"
            , "WUE_T" = "Tuber water use efficiency (g L^-1^)"
            )

4 Import data

fb <- gs %>% 
  range_read("fb") %>% 
  rename_with("tolower") %>%
  rename_with(~gsub("\\s+|\\.", "_", .)) %>% 
  mutate(treat = case_when(treat %in% "wellwater" ~ "WW"
                           , treat %in% "drydown" ~ "WD")) %>% 
  dplyr::select(block, treat, genotype,
         spad_29 = spad_29dap, 
         spad_59 = spad_59dap, 
         spad_76 = spad_76dap, 
         spad_83 = spad_83dap,
         hgt = hgt_86dap,
         rwc = rwc_84dap,
         lop = op_84dap,
         ldw = leafdw,
         sdw = stemdw,
         rdw = rootdw,
         tdw = tubdw,
         ntub,
         trs = ttrns,
         lfa = la,
         rtl = rlg
         ) %>% 
  mutate(tdb = (ldw+sdw+rdw+tdw),
         hi = tdw/(ldw+sdw+rdw+tdw),
         sla = lfa/ldw,
         rcc = (spad_83/lfa),
         wue_b = (ldw+sdw+rdw+tdw)/trs,
         wue_t = tdw/trs
         ) %>% 
  mutate(gnt = recode(genotype, !!!geno), .after = genotype) %>% 
  dplyr::select(genotype, gnt, treat, block, everything()) %>% 
  mutate(across(where(is.character), as.factor)) %>% 
  rename_with(., toupper, where(is.numeric)) 

# str(fb)

fb %>% web_table()

5 Abbreviations

gs %>% 
  range_read("var") %>% 
  filter(include == "x") %>%
  dplyr::select(Variable, Abbreviation) %>%
  kable()
Variable Abbreviation
Soil Plant Analysis Development SPAD
Plant height HGT
Relative water content RWC
Leaf osmotic potential LOP
Leaf dry weight LDW
Stem dry weight SDW
Root dry weight RDW
Tuber dry weight TDW
Tuber number NTUB
Root length RTL
Total transpiration TRS
Leaf area LFA
Total dry biomass TDB
Harvest index HI
Specific leaf area SLA
Water use efficiency WUE
Tuber water use efficiency TWUE

6 Table 1

tab <- gs %>% 
  range_read("gnt") %>% 
  dplyr::select(Number, Genotypes, Adaptability, "Growing period" = GPL, "Heat tolerance" = Heat, "Dry matter (%)")
  
# tab %>% sheet_write(gs, data = .,  sheet = "tab1")

tab %>% web_table()

7 Table 2

tab <- fb %>%
  dplyr::select(!c(block, gnt, genotype)) %>% 
  tbl_summary(
    by = treat, 
    statistic = list(all_continuous() ~ "{mean} ± {sd}"),
    missing = "no") %>% 
  add_p() %>% 
  bold_labels() %>% 
  as_tibble() %>% 
  mutate(`**Characteristic**` = recode(
    `**Characteristic**`
    , "__SPAD_29__" = "Chlorophyll concentration (SPAD) at 29 dap"
    , "__SPAD_59__" = "Chlorophyll concentration (SPAD) at 59 dap"
    , "__SPAD_76__" = "Chlorophyll concentration (SPAD) at 76 dap"
    , "__SPAD_83__" = "Chlorophyll concentration (SPAD) at 83 dap"
    , "__HGT__" = "Plant height (cm)"
    , "__RWC__" = "Relative water content (%)"
    , "__LOP__" = "Leaf osmotic potential (MPa)"
    , "__LDW__" = "Leaf dry weight (g)"
    , "__SDW__" = "Stem dry weight (g)"
    , "__RDW__" = "Root dry weight (g)"
    , "__TDW__" = "Tuber dry weight (g)"
    , "__NTUB__" = "Tuber number (N°)" 
    , "__TRS__" = "Total transpiration (mL)"
    , "__LFA__" = "Leaf area (cm^2^)"
    , "__RTL__" = "Root length (cm)"
    , "__TDB__" = "Total dry biomass (g)"
    , "__HI__"  = "Harvest index (HI)" 
    , "__SLA__" = "Specific leaf area (cm^2^ g^-1^)"
    , "__RCC__" = "Relative chlorophyll content (RCC)"
    , "__WUE_B__" = "Biomass water use efficiency (g L^-1^)"
    , "__WUE_T__" = "Tuber water use efficiency (g L^-1^)"
    )) %>% 
  rename(Variable = `**Characteristic**`
         , `Water deficit` = `**WD**, N = 75`
         , `Well-Watered`= `**WW**, N = 75`
         , `p-value` = `**p-value**`) 

# tab %>% sheet_write(gs, data = .,  sheet = "tab2")

tab %>% web_table()

8 Figure 1

fts <- gs %>% 
  range_read("ftsw") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = fts, -ID, -Genotype, -Treatment)

model <- fts %>% 
  inti::mean_comparison(response = "fts"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt1 <- model$comparison %>% 
  mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(
  type = "line", color = F
  , x = "day"
  , y = "fts"
  , group = "Treatment"
  , ylab = "Fraction of transpirable soil water"
  , xlab =  "Days"
  , glab = "Treatment"
  , legend = "none"
  , ylimits = c(0, 1.08, 0.1)
  , error =  "ste"
  ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88))

trns <- gs %>% 
  range_read("trans") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = trans, -ID, -Genotype, -Treatment) %>%
  filter(day != "TOTAL") %>%
  drop_na()


model <- trns %>% 
  inti::mean_comparison(response = "trans"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt2 <- model$comparison %>% 
    mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(type = "line", color = F
            , x = "day"
            , y = "trans"
            , group = "Treatment"
            , ylab = "Transpiration (mL)"
            , xlab =  "Days"
            , glab = "Irrigation"
            , ylimits = c(0, 700, 100)
            , error = "ste"
            ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88)
                     ) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

plot <- ggdraw(xlim = c(0, 0.5), ylim = c(0,0.9)) +
  draw_plot(plt2, width = 0.5, height = 0.43 , x = 0.0, y = 0.47 ) +
  draw_plot(plt1, width = 0.496, height = 0.5, x = 0.004, y = 0.0 ) +
          draw_plot_label(
            label = c("a", "b")
            , size = c(16, 15)
            , x = c(0.465, 0.465)
            , y = c(0.79, 0.49)
            )

plot %>% 
  ggsave2(plot = .
          , "files/Fig1.pdf", width = 18, height = 15, units = "cm")

"files/Fig1.pdf" %>% include_graphics()

9 Figure 2

# tdw

model <- fb %>% 
  lm(TDW ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: TDW
##            Df  Sum Sq Mean Sq  F value                Pr(>F)    
## block       4  1797.3   449.3   5.0193             0.0009244 ***
## gnt        14 20614.5  1472.5  16.4490 < 0.00000000000000022 ***
## treat       1  9345.0  9345.0 104.3939 < 0.00000000000000022 ***
## gnt:treat  14  2126.7   151.9   1.6969             0.0657444 .  
## Residuals 113 10115.4    89.5                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcTDW <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(TDW, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

tdw <- mc %>% 
  plot_smr(type = "line"
          , color = F
          , x = "gnt"
          , xlab = ""
          , y = "emmean"
          , ylab = "Tuber dry weight (g)"
          , group = "treat"
          , glab = "Treatment"
          , legend = "none"
          , error = "SE"
          , ylimits = c(0, 100, 20)
          , sig = ".group"
          )

# rcc

model <- fb %>% 
  lm(RCC ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: RCC
##           Df    Sum Sq   Mean Sq  F value                Pr(>F)    
## block      4 0.0001547 0.0000387   0.9411              0.443576    
## gnt       14 0.0126889 0.0009063  22.0533 < 0.00000000000000022 ***
## treat      1 0.0059808 0.0059808 145.5241 < 0.00000000000000022 ***
## gnt:treat 14 0.0013211 0.0000944   2.2961              0.009104 ** 
## Residuals 98 0.0040276 0.0000411                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcRCC <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(RCC, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

rcc <- mc %>% 
  plot_smr(type = "line"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Relative chlorophyll content"
        , ylimits = c(0, 0.1, 0.02)
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , error = "SE"
        , sig = ".group"
        ) 

# hi

model <- fb %>% 
  lm(HI ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: HI
##            Df  Sum Sq  Mean Sq F value                Pr(>F)    
## block       4 0.11840 0.029599  8.9657         0.00000251877 ***
## gnt        14 2.65590 0.189707 57.4634 < 0.00000000000000022 ***
## treat       1 0.12163 0.121629 36.8421         0.00000001753 ***
## gnt:treat  14 0.07233 0.005167  1.5650                0.1001    
## Residuals 113 0.37305 0.003301                                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcHI <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(HI, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

hi <- mc %>% 
  plot_smr(type = "line"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Harvest Index"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 1, 0.2)
        , error = "SE"
        , sig = ".group"
        ) 

# wue_t

model <- fb %>% 
  lm(WUE_T ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: WUE_T
##            Df Sum Sq Mean Sq F value                Pr(>F)    
## block       4  28.55  7.1382 10.4548          0.0000003055 ***
## gnt        14 400.34 28.5959 41.8820 < 0.00000000000000022 ***
## treat       1   2.02  2.0164  2.9533               0.08844 .  
## gnt:treat  14  14.80  1.0572  1.5485               0.10534    
## Residuals 113  77.15  0.6828                                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcWUE_T <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(WUE_T, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

wue_t <- mc %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Tuber water use efficiency (g/L)"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 10, 2)
        , error = "SE"
        , sig = ".group"
        ) 

# spad_29

model <- fb %>% 
  lm(SPAD_29 ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: SPAD_29
##            Df  Sum Sq Mean Sq F value              Pr(>F)    
## block       4   26.27   6.567  0.7227              0.5781    
## gnt        14 2396.29 171.163 18.8364 <0.0000000000000002 ***
## treat       1   12.87  12.866  1.4159              0.2365    
## gnt:treat  14  146.78  10.484  1.1538              0.3202    
## Residuals 116 1054.07   9.087                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcSPAD_29 <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(SPAD_29, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

spad_29 <- mc %>% 
    mutate(treat  = case_when(
    treat  == "WD" ~ "Water Deficit (WD)"
    , treat  == "WW" ~ "Well Watered (WW)"
    )) %>%
  plot_smr(type = "line"
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "SPAD at 29 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend = "top"
        , ylimits = c(30, 70, 5)
        , error = "SE"
        , color = F
        , sig = ".group"
        )

# spad_83

model <- fb %>% 
  lm(SPAD_83 ~ block + gnt*treat, data = .)
  
anova(model)
## Analysis of Variance Table
## 
## Response: SPAD_83
##            Df  Sum Sq Mean Sq F value              Pr(>F)    
## block       4   87.52   21.88  1.5528             0.19163    
## gnt        14 1951.83  139.42  9.8947 0.00000000000003419 ***
## treat       1  722.05  722.05 51.2451 0.00000000007974715 ***
## gnt:treat  14  436.61   31.19  2.2134             0.01098 *  
## Residuals 116 1634.44   14.09                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcSPAD_83 <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(SPAD_83, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

spad_83 <- mc %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "SPAD at 83 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend= "none"
        , ylimits = c(30, 70, 5)
        , error = "SE"
        , sig = ".group"
        ) 

plot <- ggdraw(xlim = c(0.0, 1), ylim = c(0.0, 1)) +
  draw_plot(spad_29 , width = 0.49, height = 0.35, x = 0.0, y = 0.6) +
  draw_plot(spad_83, width = 0.49, height = 0.3, x = 0.5, y = 0.6) +
  draw_plot(rcc,  width = 0.49, height = 0.3, x = 0.0, y = 0.3) +
  draw_plot(tdw, width = 0.49, height = 0.3, x = 0.5, y = 0.3) +
  draw_plot(hi , width = 0.49, height = 0.3, x = 0.0, y = 0.0) +
  draw_plot(wue_t, width = 0.49, height = 0.3, x = 0.5, y = 0.0) +
          draw_plot_label(
            label = c("a", "b", "c", "d", "e", "f"),
            x = c(0.06, 0.56, 0.06, 0.56, 0.06, 0.56),
            y = c(0.89, 0.89, 0.59, 0.59, 0.29, 0.29))

plot %>% 
  ggsave(plot = .
         , "files/Fig2.pdf"
         , units = "cm", width = 28, height = 25)

"files/Fig2.pdf" %>% include_graphics()



geno4treat <- merge(mcSPAD_29, mcSPAD_83) %>% 
  merge(., mcRCC) %>% 
  merge(., mcTDW) %>%
  merge(., mcHI) %>%
  merge(., mcWUE_T) %>% 
  arrange(desc(treat)) %>% 
  rename(Genotype = gnt, Treatment = treat)
  
# geno4treat %>% sheet_write(gs, data = .,  sheet = "tabS1")

10 Rename variables in PCA


vars <- names(fb)

model <- 5:length(vars) %>% map(function(x) {
  
    fb %>% 
      H2cal(trait = vars[x]
            , gen.name = "genotype"
            , rep.n = 4
            , fixed.model = "0 + (1|block) + genotype"
            , random.model = "1 + (1|block) + (1|genotype)"
            , emmeans = T
            , plot_diag = T
            , outliers.rm = T
            )
})


tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
      }) %>% 
  bind_rows()

varnames <- tabsmr %>% 
  dplyr::select(trait, repeatability) %>% 
  mutate(name = paste0("(",round(repeatability, 2), ") ", trait)) %>% 
  dplyr::select(!repeatability) %>% 
  deframe()

11 Figure 3

fig <- magick::image_read("files/potatos.png") %>% 
  grid::rasterGrob() %>%
  ggsave("files/Fig3.pdf", plot =  ., width = 13, height = 13)

"files/Fig3.pdf" %>% include_graphics()

12 Figure 4

mvdata <- fb %>% 
  group_by(genotype, treat) %>% 
  summarise(across(where(is.numeric), ~mean(.x, na.rm = T))) %>% 
  ungroup() %>% 
  unite(trait, c(genotype, treat), sep = "-") 

plot <- tempfile(fileext = ".pdf")
pdf(file = plot, width = 20, height = 20)
mvdata %>% 
  dplyr::select(where(is.numeric)) %>% 
  dplyr::select(!ends_with(c("29", "59", "76"))) %>% 
  pairs.panels(x = .
               , hist.col="red"
               , pch = 21
               , stars = TRUE
               , scale = FALSE
               , lm = TRUE
               ) 
graphics.off()

fig <- magick::image_read_pdf(plot) %>% 
  grid::rasterGrob() %>%
  ggsave2("files/Fig4.pdf", plot =  ., width = 20, height = 20, units = "cm")

"files/Fig4.pdf" %>% include_graphics()

13 Figure 5

mv <- mvdata %>% 
  pivot_longer(!trait) %>% 
  mutate(name = recode(name, !!!varnames)) %>% 
  pivot_wider() %>% 
  column_to_rownames("trait") %>% 
  PCA(scale.unit = T, graph = F)

hcpc <- HCPC(mv, graph = F)

outfile <- tempfile(fileext = ".pdf")
pdf(outfile, width = 9, height = 9)
plot.HCPC(x = hcpc
          , choice = "map"
              , legend = list(x = "topright"
                              , cex = 0.6
                              , inset = 0.001
                              , box.lty=0
                              )
              , draw.tree = F
          )
graphics.off()

tree <- magick::image_read_pdf(outfile) %>% 
  grid::rasterGrob()

pca <- mv %>% 
  plot.PCA(choix = "var", autoLab = "auto", cex = 0.6)

plot <- plot_grid(pca, tree
                  , nrow = 1, labels = "auto")
plot %>% 
  ggsave2(plot = .
          , "files/Fig5.pdf"
          , units = "cm", width = 30, height = 15)

"files/Fig5.pdf" %>% include_graphics()

14 Supplementary information

15 Supplementary Table 1

h2plot <- tabsmr %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  dplyr::select(trait, mean, std, min, max, V.g, V.e, repeatability) %>% 
  mutate(across(where(is.numeric), ~ round(.x, digits = 2))) %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  rename(Trait = trait, H2 = repeatability)
  
h2plot  %>% web_table()

# h2plot %>% sheet_write(gs, data = .,  sheet = "tabS1")

16 Supplementary Figure 1

var <- get_pca_var(mv)

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=3.8*ppi, height=10*ppi, res=ppi)
corrplot(var$cor, 
         method="number",
         tl.col="black", 
         tl.srt=45,)
graphics.off()

pt1 <- png::readPNG(tmp) %>%
  grid::rasterGrob(interpolate = TRUE)

pt2 <- fviz_eig(mv, 
                addlabels=TRUE,
                hjust = 0.05,
                barfill="white",
                barcolor ="darkblue",
                linecolor ="red") + 
  ylim(0, 50) + 
  labs(
    title = "PCA - percentage of explained variances",
    y = "Variance (%)") +
  theme_minimal()

pt3 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 1, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 1 - variables contribution") 

pt4 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 2, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 2 - variables contribution") 

plot <- ggdraw(xlim = c(0.0, 1.0), ylim = c(0, 1.0))+
  draw_plot(pt1,  width = 0.4, height = 0.99, x = 0.6, y = 0.0) +  
  draw_plot(pt2,  width = 0.6, height = 0.34, x = 0.03, y = 0.66) +
  draw_plot(pt3, width = 0.6, height = 0.34, x = 0.03, y = 0.33) + 
  draw_plot(pt4, width = 0.6, height = 0.34, x = 0.03, y = 0.0) +
        draw_plot_label(
    label = c("a", "b", "c", "d"),
    x = c(0.005, 0.005, 0.005, 0.65),
    y = c(0.999, 0.67, 0.34, 0.999))

ggsave2(plot = plot, "files/FigS1.pdf", height = 25, width = 30, units = "cm")

"files/FigS1.pdf" %>% include_graphics()

17 Convert figures to pdf to png

figs <- "files/"  %>%  list.files(pattern = "pdf", full.names = T)

convert <- 1:length(figs) %>% map(\(x) {
  figname <- figs[x] %>% basename() %>% gsub("pdf", "png", .)
  
  figs[x] %>% 
    image_read_pdf() %>% 
    image_write(file.path("files", figname))
})